home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / calc / pell.cal < prev    next >
Text File  |  1995-07-17  |  1KB  |  75 lines
  1. /*
  2.  * Copyright (c) 1993 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1.
  7.  * Type the solution to pells equation for a particular D.
  8.  */
  9.  
  10. define pell(D)
  11. {
  12.     local X, Y;
  13.  
  14.     X = pellx(D);
  15.     if (isnull(X)) {
  16.         print "D=":D:" is square";
  17.         return;
  18.     }
  19.     Y = isqrt((X^2 - 1) / D);
  20.     print X : "^2 - " : D : "*" : Y : "^2 = " : X^2 - D*Y^2;
  21. }
  22.  
  23.  
  24. /*
  25.  * Function to solve Pell's equation
  26.  * Returns the solution X to:
  27.  *    X^2 - D * Y^2 = 1
  28.  */
  29. define pellx(D)
  30. {
  31.     local R, Rp, U, Up, V, Vp, A, T, Q1, Q2, n;
  32.     local mat ans[2,2];
  33.     local mat tmp[2,2];
  34.  
  35.     R = isqrt(D);
  36.     Vp = D - R^2;
  37.     if (Vp == 0)
  38.         return;
  39.     Rp = R + R;
  40.     U = Rp;
  41.     Up = U;
  42.     V = 1;
  43.     A = 0;
  44.     n = 0;
  45.     ans[0,0] = 1;
  46.     ans[1,1] = 1;
  47.     tmp[0,1] = 1;
  48.     tmp[1,0] = 1;
  49.     do {
  50.         T = V;
  51.         V = A * (Up - U) + Vp;
  52.         Vp = T;
  53.         A = U // V;
  54.         Up = U;
  55.         U = Rp - U % V;
  56.         tmp[0,0] = A;
  57.         ans *= tmp;
  58.         n++;
  59.     } while (A != Rp);
  60.     Q2 = ans[[1]];
  61.     Q1 = isqrt(Q2^2 * D + 1);
  62.     if (isodd(n)) {
  63.         T = Q1^2 + D * Q2^2;
  64.         Q2 = Q1 * Q2 * 2;
  65.         Q1 = T;
  66.     }
  67.     return Q1;
  68. }
  69.  
  70. global lib_debug;
  71. if (lib_debug >= 0) {
  72.     print "pell(D) defined";
  73.     print "pellx(D) defined";
  74. }
  75.